I use the spaero::get_stats() function to calculate EWS over a moving window for the measles case incidence data from four cities in Niger. I use a bandwidth of 35 weeks, resulting in a 70 week window. All EWS are calculated with backward_only = TRUE so that EWS are only calculated based on data in the past.
# Load packages -----------------------------------------------------------
library(tidyverse)
library(spaero)
# Load data ---------------------------------------------------------------
fname <- "../../data/clean-data/weekly-measles-incidence-niger-cities-clean.RDS"
measles_data <- readRDS(fname) %>%
filter(year > 1994) # drop first NA year, only used for modeling
# Calculate EWS -----------------------------------------------------------
all_stats <- tibble() # empty tibble for storage
for(do_region in unique(measles_data$region)){
cases <- measles_data %>%
filter(region == do_region) %>%
pull(cases)
city_stats <- spaero::get_stats(
x = cases,
center_trend = "local_constant",
center_kernel = "uniform",
center_bandwidth = 35,
stat_trend = "local_constant",
stat_kernel = "uniform",
stat_bandwidth = 35,
lag = 1,
backward_only = TRUE
)$stats
city_stats_tb <- as_tibble(city_stats) %>%
mutate(
time_iter = 1:n(),
date = unique(measles_data$date)
) %>%
gather(key = ews, value = value, -time_iter, -date) %>%
mutate(region = do_region)
all_stats <- bind_rows(all_stats, city_stats_tb)
}
Define a function to plot the EWS and the observations.
plot_it <- function(dates, cases, ews, lab, title){
par(mar = c(5,5,2,5))
plot(dates, cases, type = "l", main = title)
par(new = T)
plot(dates, ews, type = "l", col = "red", axes=F, xlab=NA, ylab=NA)
axis(side = 4)
mtext(side = 4, line = 3, lab, cex = 0.7)
}
Agadez
do_region <- "Agadez (City)"
dates <- unique(measles_data$date)
par(mfrow = c(4,3))
for(do_ews in unique(all_stats$ews)){
cases <- filter(measles_data, region == do_region) %>% pull(cases)
ews <- filter(all_stats, region == do_region & ews == do_ews) %>% pull(value)
plot_it(dates, cases, ews, lab = "ews", title = do_ews)
}

Maradi
do_region <- "Maradi (City)"
dates <- unique(measles_data$date)
par(mfrow = c(4,3))
for(do_ews in unique(all_stats$ews)){
cases <- filter(measles_data, region == do_region) %>% pull(cases)
ews <- filter(all_stats, region == do_region & ews == do_ews) %>% pull(value)
plot_it(dates, cases, ews, lab = "ews", title = do_ews)
}

Niamey
do_region <- "Niamey (City)"
dates <- unique(measles_data$date)
par(mfrow = c(4,3))
for(do_ews in unique(all_stats$ews)){
cases <- filter(measles_data, region == do_region) %>% pull(cases)
ews <- filter(all_stats, region == do_region & ews == do_ews) %>% pull(value)
plot_it(dates, cases, ews, lab = "ews", title = do_ews)
}

Zinder
do_region <- "Zinder (City)"
dates <- unique(measles_data$date)
par(mfrow = c(4,3))
for(do_ews in unique(all_stats$ews)){
cases <- filter(measles_data, region == do_region) %>% pull(cases)
ews <- filter(all_stats, region == do_region & ews == do_ews) %>% pull(value)
plot_it(dates, cases, ews, lab = "ews", title = do_ews)
}

Post-outbreak EWS
As the plots above show, most EWS get swamped out by the “memory” of outbreaks. To try and get around this, I am going to look at EWS after resetting them following outbreaks. I don’t think this is too contrived because real-world early detection systems would have to do the same thing.
How to ID outbreaks
Following Toby’s suggestion, I am going to look at the distribution of outbreak sizes (i.e., cumulative cases in each year).
pop_data <- tibble(
region = c("Agadez (City)", "Maradi (City)", "Niamey (City)", "Zinder (City)"),
pop = c(118224, 267249, 1027000, 322935)
)
cusum_data <- measles_data %>%
group_by(region, year) %>%
summarise(cases = sum(cases)) %>%
left_join(pop_data)
Joining, by = "region"
ggplot(cusum_data, aes(x = cases/pop)) +
geom_histogram(aes(fill = region))

cusum_cases <- sort(cusum_data$cases / cusum_data$pop)
fit1 <- MASS::fitdistr(cusum_cases, "exponential")
x <- cusum_cases
mixexplik <- function(p,lambda1,lambda2) {
z <- p*dexp(x,lambda1) + (1-p)*dexp(x,lambda2)
return(-sum(log(z)))
}
mle_fit <- stats4::mle(minuslogl=mixexplik,
start= list(p=0.2,
lambda1 = as.numeric(fit1$estimate),
lambda2 = as.numeric(fit1$estimate)),
method="Nelder-Mead")
NaNs producedNaNs producedNaNs produced
z1 <- rexp(1000, rate = mle_fit@coef["lambda1"])
z2 <- rexp(1000, rate = mle_fit@coef["lambda2"])
zall <- mle_fit@coef["p"]*z1 + (1-mle_fit@coef["p"])*z2
plot(density(zall))
lines(density(z1), col = "red")
lines(density(z2), col = "blue")

hist(x, breaks = 25, freq = FALSE, main = NULL, ylim = c(0,1300),
xlab = "incidence (reports/person)", col = "grey")
lines(density(z1), col = "red", lwd = 2)
lines(density(z2), col = "blue", lwd = 2)
legend(0.0025, 800, legend=c("small outbreaks", "large outbreaks"),
col=c("red", "blue"), lty = 1, cex = 0.8, lwd = 2)

LS0tCnRpdGxlOiAiRVdTIGluIG1lYXNsZXMgZGF0YSIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKSSB1c2UgdGhlIGBzcGFlcm86OmdldF9zdGF0cygpYCBmdW5jdGlvbiB0byBjYWxjdWxhdGUgRVdTIG92ZXIgYSBtb3Zpbmcgd2luZG93IGZvciB0aGUgbWVhc2xlcyBjYXNlIGluY2lkZW5jZSBkYXRhIGZyb20gZm91ciBjaXRpZXMgaW4gTmlnZXIuCkkgdXNlIGEgYmFuZHdpZHRoIG9mIDM1IHdlZWtzLCByZXN1bHRpbmcgaW4gYSBgciAzNSoyYCB3ZWVrIHdpbmRvdy4KQWxsIEVXUyBhcmUgY2FsY3VsYXRlZCB3aXRoIGBiYWNrd2FyZF9vbmx5ID0gVFJVRWAgc28gdGhhdCBFV1MgYXJlIG9ubHkgY2FsY3VsYXRlZCBiYXNlZCBvbiBkYXRhIGluIHRoZSBwYXN0LgoKYGBge3IgY2FsYy1ld3N9CiMgTG9hZCBwYWNrYWdlcyAtLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLQoKbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkoc3BhZXJvKQoKCiMgTG9hZCBkYXRhIC0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLQoKZm5hbWUgPC0gIi4uLy4uL2RhdGEvY2xlYW4tZGF0YS93ZWVrbHktbWVhc2xlcy1pbmNpZGVuY2UtbmlnZXItY2l0aWVzLWNsZWFuLlJEUyIKbWVhc2xlc19kYXRhIDwtIHJlYWRSRFMoZm5hbWUpICU+JQogIGZpbHRlcih5ZWFyID4gMTk5NCkgICMgZHJvcCBmaXJzdCBOQSB5ZWFyLCBvbmx5IHVzZWQgZm9yIG1vZGVsaW5nCgoKIyBDYWxjdWxhdGUgRVdTIC0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tCgphbGxfc3RhdHMgPC0gdGliYmxlKCkgICMgZW1wdHkgdGliYmxlIGZvciBzdG9yYWdlCgpmb3IoZG9fcmVnaW9uIGluIHVuaXF1ZShtZWFzbGVzX2RhdGEkcmVnaW9uKSl7CiAgCiAgY2FzZXMgPC0gbWVhc2xlc19kYXRhICU+JQogICAgZmlsdGVyKHJlZ2lvbiA9PSBkb19yZWdpb24pICU+JQogICAgcHVsbChjYXNlcykKICAKICBjaXR5X3N0YXRzIDwtIHNwYWVybzo6Z2V0X3N0YXRzKAogICAgeCA9IGNhc2VzLAogICAgY2VudGVyX3RyZW5kID0gImxvY2FsX2NvbnN0YW50IiwgCiAgICBjZW50ZXJfa2VybmVsID0gInVuaWZvcm0iLCAKICAgIGNlbnRlcl9iYW5kd2lkdGggPSAzNSwgCiAgICBzdGF0X3RyZW5kID0gImxvY2FsX2NvbnN0YW50IiwgCiAgICBzdGF0X2tlcm5lbCA9ICJ1bmlmb3JtIiwgCiAgICBzdGF0X2JhbmR3aWR0aCA9IDM1LCAKICAgIGxhZyA9IDEsIAogICAgYmFja3dhcmRfb25seSA9IFRSVUUKICApJHN0YXRzCiAgCiAgY2l0eV9zdGF0c190YiA8LSBhc190aWJibGUoY2l0eV9zdGF0cykgJT4lCiAgICBtdXRhdGUoCiAgICAgIHRpbWVfaXRlciA9IDE6bigpLAogICAgICBkYXRlID0gdW5pcXVlKG1lYXNsZXNfZGF0YSRkYXRlKQogICAgKSAlPiUKICAgIGdhdGhlcihrZXkgPSBld3MsIHZhbHVlID0gdmFsdWUsIC10aW1lX2l0ZXIsIC1kYXRlKSAlPiUKICAgIG11dGF0ZShyZWdpb24gPSBkb19yZWdpb24pCiAgCiAgYWxsX3N0YXRzIDwtIGJpbmRfcm93cyhhbGxfc3RhdHMsIGNpdHlfc3RhdHNfdGIpCn0KCmBgYAoKRGVmaW5lIGEgZnVuY3Rpb24gdG8gcGxvdCB0aGUgRVdTIGFuZCB0aGUgb2JzZXJ2YXRpb25zLgoKYGBge3IgcGxvdC1ld3MtZnVuY30KcGxvdF9pdCA8LSBmdW5jdGlvbihkYXRlcywgY2FzZXMsIGV3cywgbGFiLCB0aXRsZSl7CiAgcGFyKG1hciA9IGMoNSw1LDIsNSkpCiAgcGxvdChkYXRlcywgY2FzZXMsIHR5cGUgPSAibCIsIG1haW4gPSB0aXRsZSkKICBwYXIobmV3ID0gVCkKICBwbG90KGRhdGVzLCBld3MsIHR5cGUgPSAibCIsIGNvbCA9ICJyZWQiLCBheGVzPUYsIHhsYWI9TkEsIHlsYWI9TkEpCiAgYXhpcyhzaWRlID0gNCkKICBtdGV4dChzaWRlID0gNCwgbGluZSA9IDMsIGxhYiwgY2V4ID0gMC43KQp9CmBgYAoKIyMjIEFnYWRlegoKYGBge3IgYWdhZGV6LCBmaWcuaGVpZ2h0PTh9CmRvX3JlZ2lvbiA8LSAiQWdhZGV6IChDaXR5KSIKZGF0ZXMgPC0gdW5pcXVlKG1lYXNsZXNfZGF0YSRkYXRlKSAKcGFyKG1mcm93ID0gYyg0LDMpKQpmb3IoZG9fZXdzIGluIHVuaXF1ZShhbGxfc3RhdHMkZXdzKSl7CiAgY2FzZXMgPC0gZmlsdGVyKG1lYXNsZXNfZGF0YSwgcmVnaW9uID09IGRvX3JlZ2lvbikgJT4lIHB1bGwoY2FzZXMpCiAgZXdzIDwtIGZpbHRlcihhbGxfc3RhdHMsIHJlZ2lvbiA9PSBkb19yZWdpb24gJiBld3MgPT0gZG9fZXdzKSAlPiUgcHVsbCh2YWx1ZSkKICBwbG90X2l0KGRhdGVzLCBjYXNlcywgZXdzLCBsYWIgPSAiZXdzIiwgdGl0bGUgPSBkb19ld3MpCn0KYGBgCgojIyMgTWFyYWRpCmBgYHtyIG1hcmFkaSwgZmlnLmhlaWdodD04fQpkb19yZWdpb24gPC0gIk1hcmFkaSAoQ2l0eSkiCmRhdGVzIDwtIHVuaXF1ZShtZWFzbGVzX2RhdGEkZGF0ZSkgCnBhcihtZnJvdyA9IGMoNCwzKSkKZm9yKGRvX2V3cyBpbiB1bmlxdWUoYWxsX3N0YXRzJGV3cykpewogIGNhc2VzIDwtIGZpbHRlcihtZWFzbGVzX2RhdGEsIHJlZ2lvbiA9PSBkb19yZWdpb24pICU+JSBwdWxsKGNhc2VzKQogIGV3cyA8LSBmaWx0ZXIoYWxsX3N0YXRzLCByZWdpb24gPT0gZG9fcmVnaW9uICYgZXdzID09IGRvX2V3cykgJT4lIHB1bGwodmFsdWUpCiAgcGxvdF9pdChkYXRlcywgY2FzZXMsIGV3cywgbGFiID0gImV3cyIsIHRpdGxlID0gZG9fZXdzKQp9CmBgYAoKIyMjIE5pYW1leQpgYGB7ciBuaWFtZXksIGZpZy5oZWlnaHQ9OH0KZG9fcmVnaW9uIDwtICJOaWFtZXkgKENpdHkpIgpkYXRlcyA8LSB1bmlxdWUobWVhc2xlc19kYXRhJGRhdGUpIApwYXIobWZyb3cgPSBjKDQsMykpCmZvcihkb19ld3MgaW4gdW5pcXVlKGFsbF9zdGF0cyRld3MpKXsKICBjYXNlcyA8LSBmaWx0ZXIobWVhc2xlc19kYXRhLCByZWdpb24gPT0gZG9fcmVnaW9uKSAlPiUgcHVsbChjYXNlcykKICBld3MgPC0gZmlsdGVyKGFsbF9zdGF0cywgcmVnaW9uID09IGRvX3JlZ2lvbiAmIGV3cyA9PSBkb19ld3MpICU+JSBwdWxsKHZhbHVlKQogIHBsb3RfaXQoZGF0ZXMsIGNhc2VzLCBld3MsIGxhYiA9ICJld3MiLCB0aXRsZSA9IGRvX2V3cykKfQpgYGAKCiMjIyBaaW5kZXIKYGBge3IgemluZGVyLCBmaWcuaGVpZ2h0PTh9CmRvX3JlZ2lvbiA8LSAiWmluZGVyIChDaXR5KSIKZGF0ZXMgPC0gdW5pcXVlKG1lYXNsZXNfZGF0YSRkYXRlKSAKcGFyKG1mcm93ID0gYyg0LDMpKQpmb3IoZG9fZXdzIGluIHVuaXF1ZShhbGxfc3RhdHMkZXdzKSl7CiAgY2FzZXMgPC0gZmlsdGVyKG1lYXNsZXNfZGF0YSwgcmVnaW9uID09IGRvX3JlZ2lvbikgJT4lIHB1bGwoY2FzZXMpCiAgZXdzIDwtIGZpbHRlcihhbGxfc3RhdHMsIHJlZ2lvbiA9PSBkb19yZWdpb24gJiBld3MgPT0gZG9fZXdzKSAlPiUgcHVsbCh2YWx1ZSkKICBwbG90X2l0KGRhdGVzLCBjYXNlcywgZXdzLCBsYWIgPSAiZXdzIiwgdGl0bGUgPSBkb19ld3MpCn0KYGBgCgojIFBvc3Qtb3V0YnJlYWsgRVdTCgpBcyB0aGUgcGxvdHMgYWJvdmUgc2hvdywgbW9zdCBFV1MgZ2V0IHN3YW1wZWQgb3V0IGJ5IHRoZSAibWVtb3J5IiBvZiBvdXRicmVha3MuClRvIHRyeSBhbmQgZ2V0IGFyb3VuZCB0aGlzLCBJIGFtIGdvaW5nIHRvIGxvb2sgYXQgRVdTIGFmdGVyIHJlc2V0dGluZyB0aGVtIGZvbGxvd2luZyBvdXRicmVha3MuCkkgZG9uJ3QgdGhpbmsgdGhpcyBpcyB0b28gY29udHJpdmVkIGJlY2F1c2UgcmVhbC13b3JsZCBlYXJseSBkZXRlY3Rpb24gc3lzdGVtcyB3b3VsZCBoYXZlIHRvIGRvIHRoZSBzYW1lIHRoaW5nLgoKIyMgSG93IHRvIElEIG91dGJyZWFrcwpGb2xsb3dpbmcgVG9ieSdzIHN1Z2dlc3Rpb24sIEkgYW0gZ29pbmcgdG8gbG9vayBhdCB0aGUgZGlzdHJpYnV0aW9uIG9mIG91dGJyZWFrIHNpemVzIChpLmUuLCBjdW11bGF0aXZlIGNhc2VzIGluIGVhY2ggeWVhcikuCgpgYGB7cn0KcG9wX2RhdGEgPC0gdGliYmxlKAogIHJlZ2lvbiA9IGMoIkFnYWRleiAoQ2l0eSkiLCAiTWFyYWRpIChDaXR5KSIsICJOaWFtZXkgKENpdHkpIiwgIlppbmRlciAoQ2l0eSkiKSwKICBwb3AgPSBjKDExODIyNCwgMjY3MjQ5LCAxMDI3MDAwLCAzMjI5MzUpCikKCmN1c3VtX2RhdGEgPC0gbWVhc2xlc19kYXRhICU+JQogIGdyb3VwX2J5KHJlZ2lvbiwgeWVhcikgJT4lCiAgc3VtbWFyaXNlKGNhc2VzID0gc3VtKGNhc2VzKSkgJT4lCiAgbGVmdF9qb2luKHBvcF9kYXRhKQoKZ2dwbG90KGN1c3VtX2RhdGEsIGFlcyh4ID0gY2FzZXMvcG9wKSkgKwogIGdlb21faGlzdG9ncmFtKGFlcyhmaWxsID0gcmVnaW9uKSkKCmN1c3VtX2Nhc2VzIDwtIHNvcnQoY3VzdW1fZGF0YSRjYXNlcyAvIGN1c3VtX2RhdGEkcG9wKQpmaXQxIDwtIE1BU1M6OmZpdGRpc3RyKGN1c3VtX2Nhc2VzLCAiZXhwb25lbnRpYWwiKSAKCnggPC0gY3VzdW1fY2FzZXMKbWl4ZXhwbGlrIDwtIGZ1bmN0aW9uKHAsbGFtYmRhMSxsYW1iZGEyKSB7CiAgeiA8LSBwKmRleHAoeCxsYW1iZGExKSArICgxLXApKmRleHAoeCxsYW1iZGEyKQogIHJldHVybigtc3VtKGxvZyh6KSkpCn0KCm1sZV9maXQgPC0gc3RhdHM0OjptbGUobWludXNsb2dsPW1peGV4cGxpaywgCiAgICAgICAgICAgICAgICAgICAgICAgc3RhcnQ9IGxpc3QocD0wLjIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbGFtYmRhMSA9IGFzLm51bWVyaWMoZml0MSRlc3RpbWF0ZSksIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGxhbWJkYTIgPSBhcy5udW1lcmljKGZpdDEkZXN0aW1hdGUpKSwgCiAgICAgICAgICAgICAgICAgICAgICAgbWV0aG9kPSJOZWxkZXItTWVhZCIpCgp6MSA8LSByZXhwKDEwMDAsIHJhdGUgPSBtbGVfZml0QGNvZWZbImxhbWJkYTEiXSkKejIgPC0gcmV4cCgxMDAwLCByYXRlID0gbWxlX2ZpdEBjb2VmWyJsYW1iZGEyIl0pCnphbGwgPC0gbWxlX2ZpdEBjb2VmWyJwIl0qejEgKyAoMS1tbGVfZml0QGNvZWZbInAiXSkqejIKcGxvdChkZW5zaXR5KHphbGwpKQpsaW5lcyhkZW5zaXR5KHoxKSwgY29sID0gInJlZCIpCmxpbmVzKGRlbnNpdHkoejIpLCBjb2wgPSAiYmx1ZSIpCgpoaXN0KHgsIGJyZWFrcyA9IDI1LCBmcmVxID0gRkFMU0UsIG1haW4gPSBOVUxMLCB5bGltID0gYygwLDEzMDApLAogICAgIHhsYWIgPSAiaW5jaWRlbmNlIChyZXBvcnRzL3BlcnNvbikiLCBjb2wgPSAiZ3JleSIpCmxpbmVzKGRlbnNpdHkoejEpLCBjb2wgPSAicmVkIiwgbHdkID0gMikKbGluZXMoZGVuc2l0eSh6MiksIGNvbCA9ICJibHVlIiwgbHdkID0gMikKbGVnZW5kKDAuMDAyNSwgODAwLCBsZWdlbmQ9Yygic21hbGwgb3V0YnJlYWtzIiwgImxhcmdlIG91dGJyZWFrcyIpLAogICAgICAgY29sPWMoInJlZCIsICJibHVlIiksIGx0eSA9IDEsIGNleCA9IDAuOCwgbHdkID0gMikKCmBgYAoKCg==